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ABSTRACT 

Numerical simulations of hot accretion flow, both hydrodynamical and magnetohydrodynamical, 
have shown that the mass accretion rate decreases with decreasing radius; consequently the density 
profile of accretion flow becomes flatter compared to the case of a constant accretion rate. This result 
has important theoretical and observational implications. However, because of technical difficulties, 
the radial dynamic range in almost all previous simulations usually spans at most two orders of 
magnitude. This small dynamical range, combined with the effects of boundary conditions, makes 
the simulation results suspectable. Especially, the radial profiles of density and inflow rate may not 
be precise enough to be used to compare with observations. In this paper we present a "two-zone" 
approach to expand the radial dynamical range from two to four orders of magnitude. We confirm 
previous results and find that from r s to 10 r s the radial profiles of accretion rate and density can be 
well described by M(r) oc r s and p oc r~ p . The values of (s, p) are (0.48, 0.65) and (0.4, 0.85), for 
viscous parameter a — 0.001 and 0.01, respectively. Or more precisely, the accretion rate is constant 
(i.e., s = 0) within ~ 10r s ; but beyond 10r s , we have s = 0.65 and 0.54 for a — 0.001 and 0.01, 
respectively. We find that the values of both s and p are similar in all numerical simulation works, 
including previous and the present ones, no matter a magnetic field is included or not and what kind 
of initial conditions are adopted. Such an apparently surprising "common" result can be explained 
by the most updated version of the adiabatic inflow-outflow model (ADIOS). The density profile 
we obtain is in good quantitative agreement with that obtained from the detailed observations and 
modeling to Sgr A* and NGC 3115. The origin and implication of such a profile will be investigated 
in a subsequent paper. 

Subject headings: accretion, accretion discs - hydrodynamics- black hole physics 



1. INTRODUCTION 

Hot accretion flow, such as advection-dominated ac- 
cretion flows (ADAFs; Ichimura 1977; Rees et al. 1982; 
Narayan & Yi 1994; 1995; Abramowicz et al. 1995; see 
Narayan, Mahadevan & Quataert 1998; Kato, Fukue & 
Mincshige 1998; and Yuan & Narayan 2013 for reviews), 
is of great interest because of its widespread applica- 
tions in low-luminosity AGNs, including the spuermas- 
sive black hole in our Galactic center, and the quiescent 
and hard states of black hole X-ray binaries (see reviews 
by Narayan 2005; Yuan 2007; Narayan & McClintock 
2008; Ho 2008; and Yuan 2011). In the early analytical 
works, the mass accretion rate is assumed to be indepen- 
dent of radius, M(r) = constant. In this case, the radial 
profile of density satisfies p(r) oc r~ 3 / 2 (e.g., Narayan & 
Yi 1994). Numerous hydrodynamical (HD) and magne- 
hydrodynamical (MHD) numerical simulations have been 
done, with most of them focusing on the dynamics of 
the accretion flow (e.g., Igumenshchev & Abramowicz 
1999, 2000; Stone, Pringle & Begelman 1999, hereafter 
SPB99; Stone & Pringle 2001; Hawley, Balbus & Stone 
2001; Machida, Matsumoto & Mineshige 2001; McKin- 
ney & Gammie 2002; Hawley & Balbus 2002; Igumen- 
shchev, Narayan & Abramowicz 2003; Pen, Matzener & 
Wong 2003; De Villiers, Hawley & Krolik 2003; Proga 
& Begelman 2003a, 2003b; Pang et al. 2011; McKinney, 
Tchekhovskoy & Blandford 2012; Narayan et al. 2012). 
The effect of strong radiation was studied by Yuan & 



Bu (2010), and most recently by Li, Ostriker & Sunyaev 
(2012). One of the most surprising - which is perhaps 
also the most important - findings of these simulations 
is that the mass accretion rate (or precisely the inflow 
rate; refer to eq. [7] in the present paper for its def- 
inition), is found to be not a constant; but rather, it 
decreases with decreasing radius. Denoting the mass ac- 
cretion rate as M(r) oc r s , numerical simulations have 
found that s ~ 0.5 — 1 (see §4.1 for a review). Con- 
sequently, the density profile flattens compared to the 
previous p(r) oc r -15 : we now have p{r) oc r~ p with 
p < 1. Such results have obtained strong observational 
supports in the case of the supermassive black hole in 
our Galactic center, Sgr A* (Yuan, Quataert & Narayan 
2003; refer to §4.2.1 for details.), the low-luminosity AGN 
NGC 3115 (Wong et al. 2011; refer to §4.2.2 for details), 
and black holes in elliptical galaxies (Di Matteo et al. 
2000; Mushotzky et al. 2000). 

In addition to the obvious theoretical interest, the ra- 
dial profiles of accretion rate and density have also im- 
portant observational applications. This is because they 
will determine the emitted spectrum and other radiative 
features of an accretion flow (e.g., Quataert & Narayan 
1999; Yuan, Quataert & Narayan 2003). For example, in 
the case of Sgr A*, the mass accretion rate at the Bondi 
radius can be determined directly by observations. Then 
depending on the radial profile of the accretion rate (or 
more exactly the density), the radiation of the accretion 
flow is completely different. The radial profile of accre- 
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tion rate also determines the evolution of black hole mass 
and spin. In many current numerical simulations, due to 
resolution difficulty, we can at most resolve the Bondi ra- 
dius and determine the Bondi accretion rate there. Then 
the evolution of mass and spin of black holes will be de- 
termined by the fraction of the Bondi accretion rate that 
finally falls onto the black hole, which is determined by 
the radial profile of accretion rate. 

It is thus important to carefully investigate the radial 
profiles of accretion rate and density. One problem of 
almost all previous simulations is that the radial dynam- 
ical range is rather small, usually at most two orders of 
magnitude. This is technically because it would be so 
time-consuming that it is almost impossible to simulate 
an accretion flow, if the dynamical range is too large. In 
addition, as is well known the simulation results are usu- 
ally not reliable close to the boundary due to the bound- 
ary condition effects. These cast some shade on the pre- 
vious simulation results, especially the exact quantitative 
radial profiles of the physical quantities. The situation is 
even worse for MHD simulations compared to HD ones. 
Because the Alfven speed is very large in the region of 
low density and strong magnetic field, MHD simulation 
is much more expensive than HD simulations. The radial 
dynamical range is thus more limited, and it is harder to 
evolve the simulation for long time which further con- 
straints the range over which a steady state is reached. 

A large dynamical range is also useful to investigate 
the following problem. Previous works (SPB99; Igu- 
menshchev & Abramowicz 1999; Yuan & Bu 2010) have 
found that the Bernoulli parameter of most outflow in 
their simulations is indeed negativeQ. One may expect 
that outflow may not be able to escape to infinity, but 
may rejoin the accretion flow at a certain distance. The 
fact that all previous simulations never found the accu- 
mulation of matter in the accretion flow could be because 
of two possible reasons. One is that the simulation time 
is not long enough for accumulation to occur. However, 
the radial velocity of outflow is roughly in the order of 
the local Keplerian velocity (Paper II). We can therefore 
expect that if the outflow would finally rejoin the accre- 
tion flow, the timescale required for accumulation would 
be shorter than the accretion timescale by a factor of 
a. Thus this possibility is unlikely. Another possible 
reason is that the radial dynamical range is too small. 
The accumulation may occur at a radius larger than the 
outer boundary of all current simulations. To examine 
this possibility, a larger dynamical range is required. 

In the present paper we simulate two-dimensional HD 
accretion flows. Especially, a "two-zone" approach will 
be adopted which helps us to overcome the technical 
problem and achieve a large dynamic range spanning four 
orders of magnitude. It is widely believed that in reality 
the magnetic stress associated with the MHD turbulence 
driven by the magnetorotational instability (MM) trans- 
fers angular momentum (Balbus & Hawley 1991; 1998). 
Therefore, we should in principle use MHD simulation. 
In the present paper we don't include magnetic field but 
instead include an anomalous shear stress to mimic the 
magnetic stress. However, as we will describe in Paper II, 

1 In our subsequent paper (Yuan, Bu & Wu 2012, hereafter Pa- 
per II) we found that depending on the initial condition of the 
simulation, the Bernoulli parameter can be negative or positive. 



in HD and MHD accretion flows, the mechanisms of pro- 
ducing the accretion rate profile are different. One may 
ask whether the radial profile of accretion rate obtained 
in the present HD simulation is the same with that ob- 
tained in more realistic MHD simulations. In this regard, 
the recent work by Begelman (2012) gives a positive an- 
swer. The comparison between our HD simulation with 
MHD simulations also confirms this point (refer to §4.1 
for details). 

The structure of the present paper is as follows. In 
§2, we introduce the details of our "two-zone" simula- 
tion approach. The results of simulations are presented 
in §3. We find that the profiles of the mass accretion rate 
and density still follows a perfect power-law form and the 
power-law index almost remains unchanged compared to 
previous results. In §4 we compare our results with pre- 
vious HD and MHD simulation works (§4.1) and obser- 
vations (§4.2). We find broadly good agreement among 
them. The last section (§5) devotes to a summary. 

2. METHOD 

2.1. Equations 

The equations describing the hydrodynamics of accre- 
tion flows are taken from SPB99: 

dp 



dt 



P S7 



p— = -Vp - pVip 



d(e/p) 
' dt 



= -pV 



V • T, 



T 2 / p. 
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We use the spherical coordinates (r, 9, 0) to solve these 
equations. In the above equations, p, p, v, and e are the 
mass density, pressure, velocity, and internal energy den- 
sity, respectively. Here d/dt = d/dt + v- V. We adopt an 
adiabatic equation of state P — (7 — l)e with 7 = 5/3. if> 
is the gravitational potential. We adopt the Paczyhsky 
k Wiita (1980) potential ip = —GM/(r - r s ), where M 
is the center black hole mass, G is the gravitational con- 
stant and r s = 2GM / c 2 . The self gravity of the accretion 
flow is neglected. 

In a real accretion flow, the Maxwell stress associated 
with the MHD turbulence driven by the magnetorota- 
tional instability transfers the angular momentum. In 
the present work, we do not include magnetic field; in- 
stead we follow SPB99 and add the final terms in equa- 
tions (2) and (3) to represent the anomalous shear stress 
and the corresponding heating. Here T is anomalous 
stress tensor. To approximate the effect of magnetic 
stress, again following SPB99, we assume that only az- 
imuthal components of the stress are non-zero because 
MRI is driven only by the shear associated with the or- 
bital dynamics: 



T , 



Or 

p sin 6 d 
r dd 



\sm6J 



(4) 



(5) 



This approximation is supported by three-dimensional 
MHD simulations (e.g., Hawley, Gammie & Balbus 1996; 
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Stone et al. f996). We emphasize this approxima- 
tion because the inclusion of other components was 
found to produce quite different results (Igumenshchev 
& Abramowicz f999). The viscosity coefficient p, = vp 
determines the magnitude of the stress and v is the 
kinematic viscosity coefficient. We adopt the form v = 
ar 1 / 2 because it corresponds to the usual "a" descrip- 
tion (SPB99). In the present paper we investigate two 
cases, i.e., a=0.001 (Model A) and a=0.01 (Model B). 
We choose units such that M = G = 1 and c = lQy/2. 

2.2. Two-zone approach 

We aim at simulations with a large dynamical range 
of four orders of magnitude in radius. Directly simulat- 
ing such a large range is technically almost impossible. 
For this purpose, we divide the whole simulation domain 
into two zones. The inner zone is from 1.35 r s to 200r s 
while the outer one from 100 r s to 40000 r s . We simulate 
these two zones separately, following the four steps de- 
scribed below. For each zone, we use the ZEUS-2D code 
(Stone & Norman 1992) in spherical geometry to solve 
the equations. The two modifications to the code are to 
include the shear stress terms and the implementation of 
the Paczyhsky & Wiita (1980) potential. We adopt the 
same non-uniform grid in both the radial and angular 
directions as in SPB99. The resolution is N r — 168 and 
N = 88. 

Step one. We first simulate the outer zone. Follow- 
ing SPB99, the initial state of our simulation is an equi- 
librium torus with a constant specific angular momen- 
tum. The readers are referred to SPB99 for the descrip- 
tion of the torus. We set the radius of the maximum 
density at lOOOOr^, the maximum density of the torus 
Pmax = 1-0, the density of the medium p = 10~ 4 , and 
pressure po = po/r. The standard outflow boundary con- 
ditions (projection of all dynamical variables) at both the 
inner and outer boundaries are adopted. 

Step two. We then simulate the inner zone. We inject 
the gas at the outer boundary. The values of density, 
specific internal energy, and velocity of the injected flow 
are taken from the steady simulation results of Step one 
at 200r s . We do not use 100r s because we hope to avoid 
the effect of the inner boundary. We adopt the outflow 
boundary condition at the inner boundary. 

Step three. Obviously, simply connecting the simula- 
tion results of the last two steps is generally not self- 
consistent. This is because, for example, the inner zone 
will obviously produce some outflow at the outer bound- 
ary and these outflow will be injected into the outer zone. 
This effect has not been taken into account in Step one. 
We therefore simulate the outer zone again, but this time 
we use the results of Step two at 100r s as the inner 
boundary condition. Again, here we do not choose 200r s 
because we hope to avoid the effect of outer boundary. In 
this way, we can capture all the outflow produced from 
the inner region and preserve all their properties such as 
their velocity and internal energy, and thus their (neg- 
ative) Bernoulli parameter. Thus, we should be able to 
observe whether these outflow will rejoin the accretion 
flow and accumulate somewhere in the outer zone. After 
the outer zone reaches the steady state, we simulate the 
inner zone again, following the method described in Step 
two. 



Step four. We plot the radial profiles of various phys- 
ical quantities along different of the inner and outer 
zones. We by-eye look at the curves in each zone to 
judge the convergence, i.e., whether the results in each 
zone obtained in two adjacent steps are consistent with 
each other. If not, we do iteration to improve. Usu- 
ally we can get satisfactory convergence after up to three 
iterations. 

However, we note that this iteration approach does not 
guarantee a complete consistency between the two zones. 
In principle, we should transfer all information of the 
fluid at each time step from one zone to another. We 
assume steady solution is reached and only transfer part 
of the information at one snapshot. However, our ap- 
proach does transfer some information between the two 
zones, and these information is perhaps the most im- 
portant for the dynamics of accretion flow. Using this 
"two-zone" approach, we can easily simulate the accre- 
tion flow with a dynamic range spanning four orders of 
magnitude in radius. Compared to the usual "one-zone" 
approach, the calculation time costed is about lO 3 /^ ~ 
several hundreds times shorter for our problem (here N 
is the number of iterations). This approach could poten- 
tially be used in other simulation problems which require 
a large dynamical range. But we should emphasize that 
the underlying assumption of this approach is that steady 
solutions exist. This is satisfied for our problem. 

3. RESULTS 

3.1. Model A: the case of a = 0.001 

Figure 1 shows the time-averaged radial distribution of 
physical quantities near the equator of the whole region 
of the accretion flow (both the inner and outer zones). 
The results are averaged over time and angle between 
6 = 84° and 6 = 96°. The solid lines are the final results 
of Step four. For comparison aim we also show by the 
dot-dashed lines the results of Step one and Step two. 
We can see that the profiles slightly become flatter after 
convergence is achieved. The density, gas pressure, ro- 
tation velocity and radial velocity can be described by a 
power law scaling with radius, 

par- 065 , Pocr- 1 ' 7 , V<t> oc r" ' 5 , v r oc r" - 55 . 

(6) 

These scalings are measured away from the innermost re- 
gion, i.e., r > 10r s . We avoid the region within 10r s be- 
cause of two reasons. Firstly, as emphasized by Narayan 
et al. (2012) and also confirmed by the radial velocity 
plot in Figure 1, close to the black hole, the radial ve- 
locity of the flow increases inward much more rapidly 
because of the strong gravity, thus it deviates from its 
scaling extrapolated from the region of r > 10r s . Sec- 
ondly, as we will see from Figure 2 below, within ~ 10r s , 
the inflow rate is a constant and there is little outflow. 
Therefore, in this sense this region is also "special" com- 
pared to the region outside ~ 10r s . We note that these 
two "inner boundary effects" do not exit if a Newtonian 
potential is adopted. 

These results are roughly consistent with SPB99, 
where they found p oc r~ ' 5 , P oc r -15 and v r oc r~ ' 5 . 
We have done some test calculations and found that the 
small discrepancy is because we adopt the Paczyhsky & 
Wiita potential while a Newtonian potential is adopted 
in SPB99. The scaling of the radial and rotational 



4 



Yuan, Wu & Bu 



>- 
oo 




oo 
oo 

Ld 



0000 



100 


0000 


10 


0000 


1 


0000 





1000 





0100 





0010 





0001 




I 00 
R 



I 000 



10000 



10.0 



> 

< 
z 
O 

I— 
< 
I— 

o 
or 




1000 10000 



1 .0000 
0.1000 

< 0.0100 
< 

01 0.0010 
0.0001 



> 



1 




1 00 1000 10000 
R 



Fig. 1. — Radial structure of the accretion flow (two zones) of Model A after iteration. All quantities have been averaged over time and 
the polar angle between 8 = 84° and 9 = 96°. The dashed line in the plot of rotational velocity denotes Keplerian rotation at the 
equator. For comparison, the results of Step one and Step two arc shown by the dot-dashed lines. 
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Fig. 2. — Time-averaged and angle integrated mass accretion rate of Model A. The solid, dashed, and dotted lines are for the inflow rate 
Min, outflow rate M ut, and net rate M aC c, respectively. For comparison, the results of inflow rate of Step one and Step two are shown by 
the dot-dashed line. 
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velocities are also consistent with the self-similar solu- 
tion of ADAF (Narayan & Yi 19940, where we have 
v r oc r _o ' 5 ,W0 oc r -0 ' 5 . The deviations of p and P from 
the self-similar solution are significant. This is because of 
the inward decrease of mass accretion rate, as we will de- 
scribe below. We note that the scaling we obtain should 
be more reliable than previous works because of our ex- 
tremely large radial dynamical range. 

In Figure 2, we plot the time-averaged and angle- 
integrated mass accretion rate of the whole region. Fol- 
lowing SPB99, the mass inflow and outflow rates, Mm 
and M out , are defined as follows, 

M in (r) = 2vrr 2 f pmin(v r ,0)smddd, (7) 
Jo 

Mout(r) = 27rr 2 / pmax(« r ,0) sin Odd. (8) 
Jo 

The net mass accretion rate is defined as, 

MaccW = Min(r) + M„ut(r). (9) 

The rates of inflow and outflow, and the net rate are de- 
noted by the solid, dashed, and dotted lines, respectively. 
Note that the results are obtained by time-average the 
integral rather than integrating the time-averages. Also 
shown in the figure by the dot-dashed line is the inflow 
rate obtained in Step one and Step two. We see that the 
final curves of inflow rate after Step jour becomes steeper 
compared to the results of Step one and Step two. The 
net rate close to the outer boundary of the outer zone 
is not constant, indicating that it requires longer time 
to get the fully steady solution there. This is because 
the viscosity coefficient is very small so the accretion 
timescale is very long. The radial profile of the inflow 
rate from ~ 2r s to 10 4 r s can be described by 

/ n. 0.48 

M fa (r) = M in (r out ) ( — ) . (10) 
Or more precisely by, 

/ x 0.65 

M in (r) = M in (r out ) — J (11) 

V '"out / 

from ~ 10r s to 10 4 r s , while it is almost constant within 
10r s . This is different from SPB99. In that work, they 
found the inflow rate keeps decreasing inward until the 
inner boundary, and the slope is steeper than ours, which 
is M in oc r 75 ("Run K" in SPB99). Again, the reason 
is because we adopt a different gravitational potential of 
the black hole. Actually, we can see from the figure that, 
beyond 100r s where the two types of potential is basi- 
cally identical, the power-law index is ~ 0.75, in good 
consistency with SPB99. The deeper potential adopted 
in our work suppresses convection to some degree and 
makes the outflow weaker. In Paper II, we will analyze 
the origin of the inward decrease of the accretion rate. 
We find that it is because of the mass loss in the out- 
flow, and the the production of outflow is because of the 

2 ADIOS has the same scaling with ADAF for the radial and 
rotational velocities. ADIOS differs from ADAF only on the scaling 
of density and accretion rate. See eqs. (14) and (16) in Narayan 
et al. (2012). 



gradient of gas pressure, i.e, the buoyant force. When 
the gravitational force is stronger, the gas pressure force 
plays a relatively minor role. In other words, in the in- 
nermost region, due to the strong gravity the accretion 
timescale is shorter than the timescale required for the 
formation of outflow thus few outflow is produced. 

If the flat density profile of p oc r~ 65 is because of the 
inward decrease of mass accretion rate, we should expect 
that the power-law index p {p oc r~ p ) should be related to 
s (M oc r s ) by p = 1.5 — s. From equations © and (ITTj) . 
we have p — 0.65 < 1.5 — s = 0.85. The small deviation 
is mainly because the density profile depends not only 
on the inflow profile, but on the sum of the inflow and 
outflow rates, i.e., M m (r) — M out (r). The radial profile 
of outflow rate is steeper than that of the inflow rate, 

M out (r) oc r - 87 (12) 

throughout the radius or 

A/out W ocr ' 8 (13) 

for the region beyond r — 100. At large radii, r > 100r s , 
the profiles of inflow and outflow rates are almost identi- 
cal thus the relation p = 1.5 — s is better satisfied, where 
0.65 and s ~ 0.75-0.8. 
The most notable result from Figures 1 & 2 is that 
the profiles of physical quantities, such as density and 
mass inflow rate, are well described by a single power- 
law form, extending from 10 4 r s to ~ 2r s . This confirms 
previous simulation results such as SPB99, although in 
those simulations the dynamical range of the whole sim- 
ulation domain usually spans only two orders of magni- 
tude, and the steady solution is only reached within one 
order of magnitude because of the effect of outer bound- 
ary condition. This result indicates that although the 
Bernoulli parameter of most of the mass outflow is neg- 
ative (SPB99 and Yuan & Bu 2010), the outflow does 
not rejoin the accretion flow and accumulate somewhere 
within 10 4 r s . Physically, this is likely because the flow 
is viscous thus there is an associated energy flux. 

3.2. Model B: the case of a = 0.01 

We simulate Model B to examine the effect of vary- 
ing the amplitude of the shear stress. The viscosity co- 
efficient in Model B is 10 times larger than Model A. 
Igumenshchev & Abramowicz (1999) and SPB99 have 
already investigated this issue. They found that when 
the stress becomes stronger, the convective instability 
becomes weaker, in the sense that the radial profile of 
inflow rate becomes flatter. But they did not present the 
qualitative results of the profiles of density and inflow 
rate when the stress is increased. In the extreme case, 
Igumenshchev & Abramowicz (1999) found that when 
a > 0.1, the flow becomes almost laminar, and a bipolar 
outflow structure very close to the rotation axis, ~ 8r s , 
is produced. 

In the present work, for the aim of comparing with 
observations (see §4 below), we try to give more qualita- 
tive results. Figures 3 is similar to Figure 1, presenting 
the radial profiles of density, pressure, rotational veloc- 
ity, and radial velocity. Figure 4 shows the radial profiles 
of inflow and outflow rates and the net rate. From Fig- 
ure 3, we see that the equatorial density, gas pressure, 
rotation velocity and radial velocity can again be well 
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Fig. 3. — Radial structure of the accretion flow (two zones) of Model B. All quantities have been averaged over time and the polar angle 
between = 84° and 9 = 96° . The dashed line in the plot of rotational velocity v$ denotes Kcplerian rotation at the equator. The density 
and pressure profiles are slightly steeper than those in Model A. 
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Fig. 4. — Time-averaged and angle integrated mass accretion rate of Model B. The solid, dashed, and dotted lines are for the mass inflow 
rate M in , outflow rate M ou t, and the net rate M a cc, respectively. 
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described by a power-law form. The slopes of the pro- 
files of rotation and radial velocity remain the same with 
Model A. The normalization of the radial velocity is 10 
times larger, as expected. Compared to Model A, the 
profiles of density and pressure are moderately steeper, 
while that of the inflow rate becomes flatter, 

p(r) oc r-°- 85 , p(r) ocr- 1 ' 85 , (14) 

Mi„(r) oc r ' 4 . (15) 

Again more precisely, the inflow rate is described by 

M in (r) oc r ' 54 (16) 

from ~ 10r s to 10 r s , while it is almost constant within 
10r s . The outflow rate is much steeper. Beyond r = 100 
it can be described by 

M out (r) oc r ' 73 . (17) 

Again we see that the relationship of p = 1.5 — s is reason- 
ably satisfied beyond 10r s , as in the case of a = 0.001. 

We have also conducted simulations with a = 0.05. 
We found that convective outflow becomes significantly 
weaker. The density profile becomes steeper while the 
accretion rate profiles becomes flatter. 

3.3. Effect of changing initial conditions 

In Model A and B, a rotating torus is adopted as the 
initial condition of our simulation. We also study the 
cases of other initial condition. One is to inject gas from 
an outer boundary, with the properties of the injected 
gas, including the temperature, rotation and radial ve- 
locity, being determined by the self-similar solution of an 
ADAF (Narayan & Yi 1994). Another model is that the 
initial condition is expanded from the one dimensional 
global solution of ADAFs. We find that the radial pro- 
files of inflow rate and density are very similar (Paper II) . 
However, this will not be the case if the angular momen- 
tum of the injected gas in the initial condition is very low. 
In that case, the profile of accretion rate will become flat- 
ter and correspondingly density profile steeper. The full 
discussion of the effect of initial condition and boundary 
condition will be presented in a subsequent paper (Bu, 
Yuan & Wu 2012). 

3.4. The new scaling law of hot accretion flow 

Narayan & Yi (1995; see also Narayan, Mahadevan 
& Quataert 1998) presented the radial scaling of many 
quantities of the hot accretion flow based on the self- 
similar solution. These solutions are very useful to esti- 
mate approximately the properties of hot accretion flow. 
But at that time, the inward decrease of the accretion 
rate has not been found and taken into account in the 
solution. For the convenience of future use, here we 
present the new scaling law. In terms of interpreting ob- 
servations, the most prominent effect of outflow is on the 
profile of density. Therefore, compared to Narayan & Yi 
(1995), only the scaling of density-related quantities such 
as the number density, magnetic field strength, pressure, 
and viscous dissipation rate, are changed. The scaling of 
other quantities are also presented for completeness. In 
the original scaling of Narayan & Yi (1995), the param- 
eter /3 = 0.5 (= p gas / (p gas +p mag ), with p mag = B 2 /8ir). 



Here we do not specify a value for ft. The new results 
are, 

v w -1.1 x 10 n ar" 1/2 cm s" 1 
ft w 2.9 x 10 4 to-V- 3/2 s- 1 
c 2 w 1.4 x lO 20 ^ 1 cm 2 s~ 2 

n e w 6.3 x lO^a^TO^moutf-ou/ 2 ( —) cm ~ 3 

\ ^*OUt / 

B « 6.5xl0 8 (l-/3) 1 / 2 a-V2 m -i/2 rh i/2 r P/2-3/4 r _i/ 2 - P /2 q 

p w 1.7 x 10 16 a" 1 m- 1 rh out r^~ t 3/2 r- 1 - p g cm" 1 s~ 2 
q+ w 5.0 x 10 21 m - 2 rh out r- 5/2 - p r^ t 3/2 erg cm" 1 s~ 2 

t cs w 24a- 1 m out r 1 - p r p mt 3/2 (18) 

Same with Narayan & Yi (1995), all quantities are writ- 
ten in scaled units, M = m M , r is radius in unit of 

r s , M = mMEdd (Afedd = 10LEdd/c 2 ), m out is the mass 
accretion rate at the outer boundary r out . 

4. COMPARISON WITH PREVIOUS WORKS AND 
OBSERVATIONS 

4.1. Comparison with previous simulations: a common 
radial density profile ? 

SPB99 found that, for the a-description (Run K in 
SPB99) as we adopte in the present work, p(r) oc r~ 0,5 . 
Igumenshchev & Abramowicz (2000) performed two- 
dimensional HD simulations, considering a larger range 
in the parameter space spanned by a and adiabatic index 
7. They found that when 7 = 5/3, for both a = 0.01 
and 0.03 the density profile is roughly the same, i.e., 
p{r) oc r~ ' 5 . This result is approximately consistent 
with SPB99. The density profile in these two works is 
flatter than our result since a Newtonian potential was 
adopted there. McKinney & Gammie (2002) obtained 
p oc r~°- 6 and v r oc r~ 2 . But this result strongly suf- 
fers from the "inner boundary effect", since a pseudo- 
Newtonian potential was adopted and the above scaling 
was measured over 1.3r s < r < 10r s . 

Many MHD numerical simulations of hot accretion flow 
have been done after the pioneer HD works of Igumen- 
shchev & Abramowicz (1999, 2000) and SPB99. For two- 
dimensional simulations, since turbulence and accretion 
die away because there is no dynamo action to main- 
tain the poloidal magnetic field, steady state can not be 
reached at large radii. This effect, combined with the in- 
ner boundary effect, severely restricts the measurement 
of the radial profiles of dynamical quantities such as den- 
sity. For three-dimensional simulation, it is very time- 
consuming to simulate a large dynamical range. Perhaps 
due to these reasons, only a few works have presented 
radial profiles of density and inflow rate. For example, 
Stone & Pringle (2001), Hawley, Balbus & Stone (2001), 
and Hawley & Balbus (2002) did the two and three- 
dimensional MHD simulation of accretion flow. They 
all found that the mass accretion rate decreases with de- 
creasing radius, but the quantitative radial profile was 
not given. 

Machida, Matsumoto & Mineshige (2001) performed 
three-dimensional MHD simulation of hot accretion 
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flows, with a toroidal initial magnetic field configuration. 
A Newtonian potential was adopted thus their results are 
not affected by the "inner boundary effect" . They found 
that the radial structure of the accretion flow is very sim- 
ilar to that obtained by HD simulations. Two fits were 
obtained. In the "early stage" of the simulation they 
found that the radial profiles of the density and radial 
velocity can be described by p oc r~ 5 and v r oc r~ 1,5 . 
At the "late stage" of their simulation, they found mod- 
erately different results, p oc r~ - 8 and v r oc r~ 13 . While 
the profile of density is very similar to our result, their 
profile of radial velocity is much steeper. This may be be- 
cause that the averaging approach adopted in their work 
is different from ours. They average the quantity over the 
whole 9 angles (eqs. [1] and [2] in Machida, Matsumoto 
& Mineshige 2001); while our average is constrained to 
84° < 6 < 96°. This approach is expected to produce a 
quite different result only to the radial velocity, since the 
radial velocity away from the equatorial plane is much 
larger than that close to 6 ~ 90° while density is on the 
opposite. Igumenshchev, Narayan & Abramowicz (2003) 
also presented three-dimensional MHD simulation. They 
adopted Paczynsky & Wiita potential and two types of 
initial magnetic field configurations. Only in the case 
of a toroidal initial configuration the density profile was 
given, which is p oc r _1 . This profile is steeper than 
that of Machida, Matsumoto & Mineshige (2001). The 
discrepancy should not because of the inner boundary 
effect, since the measurement is over 7r s < r < 100r s . 
We speculate that one reason may be that a different po- 
tential was adopted. Another reason may be associated 
with the strong fluctuation of the radial velocity in this 
radial range (see Fig. 11 in Igumenshchev, Narayan & 
Abramowicz 2003). 

In the above-mentioned MHD simulations, the mag- 
netic field is weak. Recently attention has been paid 
to accretion flows with accumulated strong magnetic 
flux, i.e., "magnetically arrested disk" or "magneti- 
cally choked accretion flows" (Narayan, Igumenshchev 
& Abramowicz 2003; Pen et al. 2003; Pang et al. 2011; 
McKinney, Tchekhovskoy & Blandford 2012). Pen et al. 
(2003) constructed three-dimensional MHD simulation 
with huge number (1400 3 ) of grid zones. They empha- 
sized that the axial and reflection symmetries usually 
adopted in many simulations are not adopted in their 
work. The input magnetic field is different from the usual 
way. It it an admixture of random field loops and large 
loops that thread the whole simulation box. The density 
profile they measured is p oc r~°' 72 over the radial range 
of 0.03r\B ~ 0.3rB, with re is the Bondi radius. The 
large distance of the inner boundary from the black hole 
implies that this result does not suffer from the "inner 
boundary effect" . This slope is well between the results 
of our Model A and Model B. 

In another three-dimensional MHD simulation by Pang 
et al. (2011), special attention was paid to the radial den- 
sity (or equivalently accretion rate) profile of the accre- 
tion flow. A Newtonian potential was adopted in their 
work. The outer boundary is quite large, almost ten 
times larger than the Bondi radius. They conducted a 
numerical survey of parameter space, namely the mag- 
nitude of the ambient magnetic field, the radial dynami- 
cal range within Bondi radius, and the resolution of the 



Bondi scale. The found that the mass accretion rate at 
the inner boundary r ln can be well described by the fol- 
lowing power-law form, 

M Bo ndi V r B J 

with the mostly favored value of k e g 1 . Here M Bondi is 
the accretion rate at Bondi radius re and r; n ~ 0.03rB- 
Again like the case of Pen et al. (2003), the rather large 
rj n ensures that the obtained profile is not affected by the 
inner boundary effect. This result is in good agreement 
with our Model B (eq. [T6] ) . 

Mckinney, Tchekhovskoy, & Blandford (2012) simu- 
lated three-dimensional general relativistic MHD accre- 
tion flow with different thickness and two types of ini- 
tial magnetic field configuration, namely poloidal and 
toroidal. Density profiles haven been obtained for 
"thick" flows, which is p oc r~ - 6 , same for both con- 
figurations. For "thin" flows, it is p oc r~ - 7 . However, 
the measurement is not relevant to us since it is over 
6 ?, 6 ^ r i$ 15r s , where the inner boundary effect is very 
strong. 

Narayan et al. (2012) most recently presented two 
long-duration GRMHD simulations of a hot accretion 
flow, in which the accumulation of magnetic flux around 
the black hole does and does not occur ("MAD" and 
"SANE" models), respectively. The steady state is 
reached up to ~ 50r s . After taking into account the 
"inner boundary effect" and u a(r) effect", they found 
v r oc r -0 ' 5 when r > 10r s . Regarding the density profile, 
they claimed that its radial scaling is such that M(r) ~ 
constant is roughly satisfied at small radii. Note that 
this is not in conflict with our results. As we have stated 
before and shown by Figs. 2&4, M(r) ps constant when 
r < 10r s . Outside ~ 10r s the radial inflow rate profiles 
for the "MAD" and "SANE" models are M(r) in oc r 0A 
and oc r ' 65 , respectively (Ramesh Narayan, private com- 
munications). The result of the "SANE" model is in good 
consistency with our Model A (eq.[TT]) while the result of 
"MAD" model is similar to our Model B (eq.pl)]). From 
the above-mentioned works we note that the inflow rate 
profile of accretion flow with accumulated magnetic flux 
seems to be somewhat flatter than that without accumu- 
lation. This may be because that when there is accumu- 
lation of magnetic flux, the "effective a" is larger, i.e., 
"MAD" and "SANE" correspond to our Model B and A, 
respectively. 

Overall, almost all current numerical simulations give 
a similar radial density profile, with p oc ^-(O- 5 - 1 ). The 
results depend weakly on the presence or absence of mag- 
netic field (HD or MHD), viscous parameter a, strength 
of the magnetic field (weak or strong), the initial con- 
figuration of the magnetic field (toroidal or poloidal or 
tangled), and the dimension of calculation (two or three 
dimension) . If this is true, an interesting question is then: 
what causes such a result? To answer this question, 
one first needs to understand what causes the inward 
decrease of the mass accretion rate. This question will 
be studied in Paper II. Here we briefly summarize the 
results. Two scenarios have been proposed to explain 
this result. One is the convection-dominated accretion 
flow (CDAF) model (Narayan et al. 2000; Quataert & 
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Gruzinov 2000a). This model proposes that both hydro 
and MHD accretion flows are convectively unstable. In 
this case the inward decrease of accretion rate is because 
that the fluid circulates in the convective eddies. In con- 
trast to the CDAF model, the adiabatic inflow-outflow 
solution (ADIOS) model (Blandford & Begelman 1999, 
2004; Begelman 2012) suggests that the inward decrease 
of accretion rate is because of mass loss in the outflow. 
Our analysis in Paper II favors the latter. Moreover, we 
find that outflow is driven by buoyant force and mag- 
netic centrifugal force in HD and MHD accretion flows, 
respectively. 

Then it is an interesting question why different phys- 
ical mechanisms cause the roughly same profiles of ac- 
cretion rate and density? Begelman (2012) gives an an- 
swer to this question. In the early version of the ADIOS 
model (Blandford & Begelman 1999, 2004), there are 
inflowing and outflowing zones with equal but opposite 
mass fluxes that vary with radius following a power-law 
form, M oc r s . The net rate, which is the sum of in- 
flow and outflow rates, is very small. They applied the 
conservation laws of hydrodynamics to the inflow zone 
and obtained relationships among the parameters such 
as the specific angular momentum, effective binding en- 
ergy, and the value of s. In this model, the value of s is 
not determined but only constrained to be in the range 
< s < 1. One advantage of the ADIOS model, as 
claimed by the authors, is that although it was based on 
HD analysis, the approach can be fully applied to MHD 
accretion flows as well. But the disadvantage is that the 
early ADIOS model is intrinsically one-dimensional, the 
outflow zone is laminar and their streamlines are self- 
consistently patched onto the inflow zone at each radius. 
Begelman (2012) recently developed this model. The 
main improvement is that now the inflow and outflow 
zones are treated on an equal footing, and the outflow 
is permitted to be as turbulent as inflow rather than a 
laminar flow. He found that the viable range of s is now 
narrowed down from < s < 1 to s » 1. The main 
reason for the regulation of the solution is an introduc- 
tion of a conserved outward flux of energy through the 
flow, which was ignored in Blandford & Begelman (1999, 
2004). This work therefore well explains why various 
simulations obtain very similar profiles of accretion rate 
and density. Because of this reason, we have confidence 
that our HD simulation results should be similar to the 
more realistic MHD simulation results, and further can 
be compared to observations. Of course, we note that the 
quantitative result of the value of s obtained in Begelman 
(2012) is moderately larger than that obtained in simu- 
lations. In Model A, when r > 100r s where the "inner 
boundary effect" is the weakest thus more appropriate to 
be compared with Begelman (2012), we find s ~ 0.75. 

4.2. Comparisons with observations 
4.2.1. SgrA* 

Because of its proximity, the supermassive black hole 
located at the center of our Galaxy is regarded as the 
unique laboratory for the study of black hole accretion 
(see Yuan 2011 for a review). Abundant of data has 
been obtained, which puts the most strict constrains to 
the theory of black hole accretion. Among them is the 
constrain on the mass accretion rate at both the Bondi 



radius and the region close to the black hole horizon. 
The high spacial resolution of Chandra can well resolve 
the Bondi radius and infer the gas density and temper- 
ature there (Baganoff et al. 2003). The values of the 
Bondi radius and the Bondi accretion rate are then well 
determined, r Bon di ~ 10 5 r s and M Bon di ~ 10~ 5 M Q yr" 1 , 
respectively. Although the value of mass accretion rate 
is obtained by the simple analytical Bondi theory, it was 
found to be in good consistency with the detailed three- 
dimensional numerical simulations focusing on the fuel- 
ing of the supermassive black hole in Sgr A* of Cuadra 
et al. (2006). In fact, the numerical simulations obtained 
in this work is M w 3 x 10~ 6 M Q yr" 1 , only a factor of 
3 lower than MBondi- This factor is likely because of the 
inclusion of angular momentum of the accretion flow in 
the simulation. 

The accretion rate at the innermost region of the ac- 
cretion flow, on the other hand, is strongly constrained 
by the radio polarization observations (e.g., Quataert & 
Gruzinov 2000b; Macquart et al. 2006; Marrone et al. 
2007). The detected high level of linear polarization at 
frequencies higher than ~ 150GHz sets an upper limit 
to the rotation measure, which then requires that the 
accretion rate at the innermost region of the accretion 
flow must be in the range between 2 x 10 _7 M Q yr -1 and 
2 x 10~ 9 M© yr -1 , depending on the assumed configura- 
tion of the magnetic field. This implies that more than 
99% of the gas captured at the Bondi radius must be 
lost and does not fall onto the black hole. Correspond- 
ingly, the density profile of the accretion flow is signifi- 
cantly flatter than the prediction of p oc r~ 3//2 . The de- 
tailed modeling to Sgr A* presented in Yuan, Quataert 
& Narayan (2003) has shown that the radial profiles of 
accretion rate and density are described by M oc r 03 and 
p oc r _1 . 

This density profile is moderately steeper than that 
obtained in our simulation. This is likely because that 
the specific angular momentum of the accretion flow at 
the Bondi radius is significantly sub-Keplerian, <~ 0.3 
times Keplerian, as the numerical simulation has indi- 
cated (Cuadra et al. 2006). In this case, our study (Bu, 
Yuan & Wu 2012) indicates that the outflow becomes 
significantly weaker, and the density profile correspond- 
ingly steeper. 

4.2.2. NGC 3115 

NGC 3115 is an low-luminosity AGN. The mass of the 
black hole in its center is M = (1 - 2) x 10 9 M . The 
source is very dim. Chandra observations only present 
an upper limit of 10 38 ergs -1 <~ 10^ 9 LEdd (Wong et al. 
2011 and references therein). Therefore the accretion 
mode in this source must be an ADAF. At a distance 
of 9.7 Mpc, NGC 3115 is the nearest > 10 9 M black 
hole. Therefore, the angular size of the accretion flow in 
this source is very large, which is very helpful for us to 
directly detect the accretion flow. Wong et al. (2011) 
recently conducted Chandra observation to this source 

and determined that the Bondi radius is about 4 — 5 . 
The most prominent result is that for the first time, they 
have resolved the accretion flow within the Bondi radius. 
They found that the temperature is rising toward the 
black hole as expected in all accretion flow models. But 
more importantly, the radial density profile of the accre- 
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tion flow within 4 is found to be well described by a 
power-law form, 

p(r) (xr- 1 - 03 ^ ^. (20) 

This result is similar to Sgr A* and is in reasonable agree- 
ment with our numerical simulations. 

5. SUMMARY 

One of the most important finding of both HD and 
MHD numerical simulations of hot accretion flow in the 
recent years is that the mass accretion rate decreases 
with decreasing radius. Correspondingly, the density 
profile becomes flatter compared to the previous pre- 
diction when the accretion rate is constant of radius. 
One main problem with previous simulation is that be- 
cause of technical difficulty the radial dynamical range 
in simulations is usually very limited, typically span- 
ning less than two orders of magnitude. The bound- 
ary condition effect makes the radial range over which 
we can reliably measure the profile of physical quantities 
even smaller. The previous simulation results are there- 
fore somewhat suspectable. In this paper we adopt a 
"two-zone" approach to simulate an axisymmetric accre- 
tion flow, extending the dynamical range to over four 
orders of magnitude, i.e, from ~ r s to 40000r s . We 
confirm previous results that the profiles of inflow rate 
and density can be well described by power-law forms. 
Within 10r s , M- m (r) ~ const.. Beyond 10r s , the power- 
law slopes are a function of the viscous parameter a. 
For a = 0.001, they are described by p(r) cx r~ 65 
and Min(r) cx r 065 . For a — 0.01, the results are 
p(r) cx r-°- S5 and M in (r) cx r - 54 . 

We also combine from literature all available numer- 
ical simulations which have presented the radial profile 
of density, both two-dimensional and three-dimensional, 
HD and MHD. We find that all these simulations give 



somewhat similar results, p cx r -(°- 5 ~ 1 ) j anc j this is also 
consistent with our results. The diversity of the power- 
law index seems to come from the differences of the value 
of a, the gravitational potential of the black hole, the 
initial condition of the simulation, and the strength of 
the magnetic field. The rough consistency among var- 
ious simulations can be explained by the most recent 
ADIOS (adiabatic inflow-outflow solution) model (Begel- 
man 2012). In this work it was found that after con- 
sidering both the inflow and outflow zones at the equal 
footing and a conserved outward energy flux, the value 
of s is well constrained to be in a narrow range. 

The radial profiles of accretion rate and density ob- 
tained by numerical simulations are in good agreement 
with observations. In the case of Sgr A*, detailed mod- 
cling to the multi- waveband spectrum find that M{r) cx 
r - 3 and p cx r _1 . In the case of NGC 3115, for the 
first time, Chandra observations resolve the accretion 
flow within Bondi radius, and find par uo -o.2i , 
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